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I. INTRODUCTION 

Recent measurements of B c meson mass and lifetime in CDF [1] and DO [2] experiments 
are the first steps in the experimental research of heavy quarkonia with open flavor. The 
obtained experimental results are in a good agreement with the theoretical predictions for 

nn 

the B c mass |3H5|: 



m cr>F = 6 2756 ± 0.0029(stat.) ± 0.0025(sys.) GeV, 
m °° = 6.3000 ± 0.0014(stat.) ± 0.0005(sys.) GeV, 



m theor = q 25 ± Q Q3 Q eV . 



as well as for the decay time 6Ml0l|: 



t cdf = o.448i°;°i(stat.) ± 0.032(sys.) ps, 
r°° = 0.475l^(stat.) ± 0.018(sys.) ps, 

r theor = Q A8 ± Q Q5 ps> 

The B c meson was observed only in two decay modes: B c — > J/ipn and B c — > J /ip+p+v^. 
The investigation of other decay modes will be possible at at LHC, where about 10 10 events 
with B c mesons per year are expected. This huge amount of events will allow to obtain the 
information on the production cross section distributions, on the decay branching fractions, 
and in some cases, on the distributions of decay products. 

The B c systems do not have strong and electromagnetic annihilation decay modes. Due 
to this reason the exited B c systems laying below B + D threshold have the decay widths by 
two order of magnitude smaller than the widths for analogous exited states of charmonium 
and bottonium. All exited B c states after a set of radiative transitions decay into the lightest 
pseudoscalar state (0~). The lifetime of this state is comparable with the lifetimes of of B 
and D mesons and essentially differ from lifetimes of other lightest quarkonia: f] c and rjb- 
This is why B c meson provides a unique possibility to investigate the both strong and weak 
interactions. 

The main B c decay modes, such as B c — > J/i/j + Iv, B r — > J/ipir and B c — > B^ p were 



theoretically studied in details (see, for example p, O, ll2|]). For all these processes it is 



assumed that the factorization approach is valid: the decay B c — > heavy hadron + W* is 
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followed by the decay of the virtual H^-boson. The transition form-factors for the processes 
B c — y heavy hadron + W* can be estimated within QCD sum rules, within quark potential 
models, or in the framework of light-cone quark model. 

nn 

In our recent articles [13|, H4J we have studied the B c decay process with several pions in 

(*) 

the final state, such as B c — y J/tp+nTr and B c — y B s +nir with 1 < n < 4. The factorization 
approach have been used in these calculations. The characteristics of the virtual M^-boson 
decay have been adopted from the experimental data on r — y v T + nn decays for n < 3 and 
from the experimental data on it mesons production in the process e + e~ —y An. 

The modules for calculation of these decay amplitude is implemented into EvtGen pro- 
gram package 15]. Therefore the detailed simulation of these decays in the kinematical 
condition of real experiments is possible now. In this paper we describe the installation and 
using the developed modules. 

This paper is organized as follows. In the next section we present a theoretical background 
for calculation of B c —y V(p) + mr decays. In sections [TTT] the structure of presented packages 
is described. The installation instruction and the comparison with the available data are 
given in section IIVI The conclusions are given in the final section. 



II. THEORETICAL BACKGROUND 



To calculate the decay amplitudes we assume that the discussed processes can be repre- 
sented as the decay B c —y heavy hadron + W* followed by the decay of the virtual VT-boson. 
An additional assumption should be made that the only one of the constituent quarks 
weakly decays, meanwhile the other quark remains the same (a typical diagram is presented 
in Fig. EE]). 

Within this approach the amplitude of the process can be written in the form 

A[B C J/i, + mr] = G^a^n'/*^ (i) 

where ew is the polarization vector of W*; 

H 3 ^ is the B c — > J/tp + W* transition vertex: 

= (J/if> |c7„ (1 - 75) b\ B c ) = V M - A,. (2) 

Vector and axial currents are equal to 

V, = (J/i/> |c 7m 6| B c ) = i^etpaqpFv (g 2 ) , (3) 
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Figure 1: Typical diagram for B c — > V{P) + nir decay. 



= (J/ip \cy^ 5 b\ B c ) 



« (? 2 ) + P, (^b c ) F + A (g 2 ) + % (e* PBc ) (q 2 ) , (4) 



where ps c and pji^ are the momenta of B c - and .//^-mesons; 
g = Pb c — Pji^ is the momentum of virtual PF-boson; 

P = Pb c + Pj/ip] 

is the polarization vector of J/ip meson; 
and F v (q 2 ), F A (g 2 ), F^(q 2 ), F A (q 2 ) and F v (q 2 ) are form-factors of B c -> J/^ + VF* 
decays. 

In the tree approximation the parameter a\ (lir) is equal to unity. Higher-order corrections 
lead to dependence of this factor on renormalization scale lir [16|. Numerical values for 
eti (lir) at different scale are calculated in . For the process B c — > J/ip + n7r the value of 
lir has been chosen to be equal to the mass value of the decayed 6-quark: 



ai(m b ) = 1.4 



The form-factors were calculated within different nonperturbative approaches: QCD sum 
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B c — > J/ip + nir 




Fi(0) 


Cl 


C2 


A 


5.9 


0.049 


0.0015 


A + 


-0.074 


0.049 


0.0015 


V 


0.11 


0.049 


0.0015 


B c ^ B* + nir 




*l(0) 


Cl 


C2 


A 


8.1 


0.30 


0.069 


A + 


0.15 


0.30 


0.069 


V 


1.08 


0.30 


0.069 


B c —J- B s + nir 




/i(0) 


Cl 


C2 


u 


1.3 


0.30 


0.069 



Table I: Form-factor parameters for different SR form-factor sets. 
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potential quark models [9|, and Light-Front quark models |l7l Il8|. In our article 



we use exponential parametrization of these form-factors: 

F^q 2 ) = F i ^)e W {c l q 2 + c 2 q i }. 



(5) 



The reason for this choice is that in the framework of potential quark models the form- 
factors of B c — > J/ip transitions are expressed through the interception integral of initial 
and final meson wave functions. For heavy quarkonia these wave functions are usually 
parametrized in exponential form, so presented above expression for form-factors is natu- 
ral. However, it can be shown, that the numerical results do not depend strongly on the 
parametrization form. In the Tab. [J we give numerical values of form-factor parameters 
Fi(0), Ci )2 for different form-factor sets. 

The decay B c — > B* + nir is described by analogy with the decay B c — > J/ip + rm. The 
same formula for the amplitude is used with the other parameter set (see Tab. [T|). 
The amplitude for the decay B c — )■ B s + nix can be written in more simple form 

G F V CS 



A[B C B s + nir} 



V2 



'-a 1 (fi R )'H 



(6) 
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where 

= (B s \c (1 - 75) b\ B c ) = f + (q 2 ) Pfl + f4q 2 )q»- (7) 

The detailed information about W* decay is not needed to obtain the integrated branching 
fractions, as well as the branching fraction distributions on q 2 . Let us consider the decay 
B c — > J/ipW* J/ipnir as an example: 

dT (B c -> J/^TZ) = ^^a\WU w e™\™d§ (B c J/^nrr) , (8) 
where Lorentz-invariant phase space is defined according to 

d$ (Q^ Pl ... Pn ) = (2tt)^ (Q - J» II 2^3- ( 9 ) 
Using the following recurrent expression for the phase space 

da 2 

rf$ (B c -»■ J/?/W* J/^vr) = — d$ (B c -»■ J/^W) rf$ (W* -»■ mr) (10) 

one can perform the integration over phase space of the final state nir: 

j-Jd$ (W* nvr) ejfe^* = (ff^ - <?V) Pr (g 2 ) + W l (g 2 ) , (11) 

where spectral functions pt,l (g 2 ) are universal and can be determined from theoretical and 
experimental analysis of some other processes, for example r — > v T rm decay or electron- 
positron annihilation e + e~ — > mr. (See fisl. [h! I19I] for details.) It is worth to mention, that 
due to the vector current conservation and the partial axial current conservation spectral 
function pi for n > 2 is negligible in almost whole kinematical region, so it can be neglected 
in the estimations for n > 2. For the purposes of our article, however, a more detailed 
description is required for the multipion final state. 

In the framework of resonance model the decays of W-boson can be described in terms 
of virtual p- and ai-mesons exchange (see typical diagrams presented in Fig. |2]). 

If there is one tt- meson in the final state (see Fig. Wfr), the vertex of W* — > it transition 
can be written in the form 

(n + \(du) V - A \W) = Ukp, (12) 

where k is the momentum of tt- meson and f n ~ 140 MeV is its coupling constant. In 
accordance with this interaction vertex the effective polarization vector ej^ in this case has 
the form 

e)f = UKI<- (13) 
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Figure 2: Typical diagrams for W* — > nir decay in resonance approximation. 

Note, that this vertex violates the axial current. 

The decays of virtual PF-boson into multipion final states are described within resonance 
model in terms of virtual p- and ai-mesons exchange (see typical diagrams presented in 
Fig.Gb,c). 

The W* — > 7r + 7r° decay is saturated mainly by contributions of virtual p- and p'-mesons 
(see Fig. [2h>). The corresponding effective polarization vector can be written as 

ef = F p (q 2 )(k 1 -k 2 ) p , (14) 

where k\ j2 are 7r-mesons momenta and F p (q 2 ) is the p-meson form-factor (see 20] ). The 
difference in 7r°- and 7r + -meson masses is neglected, thus the virtual W boson in this decay 
has a transverse polarization. It should be noted that the width of the p meson must be 
taken into account. 

The W* — > 37r-transition is saturated mainly by W* — > a\ — > pn — > 3tt decay chain. 
Following 20| one can write the effective polarization vertex in this case as 

# = -iY^F a {q 2 ) {B P (s 2 )V lp + B p { Sl )V 2 ,} , (15) 

where 



V\,2n = k Xt 2n - k 3p - q p — (16) 

and 

si,2 = (q~ K 2 f. (17) 



Parametrization of B p (s) function is presented in [20] . It can be clearly seen, that if one 
neglects the difference between charged and neutral 7r-meson masses the above expression 
in transverse and the axial current is conserved. In EvtGen package this transition was 
realized already in TAUHADNU model. 
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The simulation of W* — > 47r-transition will be possible within the next version of our 
package. 

III. OVERVIEW OF THE SOFTWARE STRUCTURE 

For the Monte-Carlo simulation of the discussed decays the generator EvtGen has been 
chosen. This generator is widely used in high energy physics. For example, EvtGen is a 
part of the LHCb software environment GAUSS. 

The characteristic property of EvtGen is that, in contrast to many other generators 
(Pythia, for example), the polarizations of the initial and final particles are taken into ac- 
count in modeling. This is possible because EvtGen deals with the amplitude of the process 
under consideration, while other generators work with averaged over helicities squared ma- 
trix elements. As a result, one can easily study different spin asymmetries and other similar 
quantities. 

EvtGen-generator is written completely on C++. In order to add a new decay to it one 
should just create a new class, that describes this decay. The base class is EvtDecayAmp, 
where prototypes of all necessary for the model functions are given. The most significant of 
these functions are init(), initProbMax() and decay() (see Fig. [3]). 

The method init() performs the initialization of the decay model and reads its parame- 
ters. The necessary parameters are stored in the so-called .dec-file (the exact name of this 
file can be determined by the user ). T he detailed description of the .dec-file format can be 



151 ] . Briefly, the structure of this file is 



found in EvtGen documentation 
BeginDecay (decP article) 

(Br) (outi) . . . (out n ) (modelName) (pari) . . . (par m ); 

EndDecay 

where 

(decName) is the name of the initial particle, 

(Br) is the branching fraction of the decay, 

(outi) . . . (out n ) are the names of the final particles, 
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EvtDecayBase 



I 

EvtDecayAmp 



I 



EvtBcNPi 


EvtBcNPi. cpp 
EvtBcNPi. hh 


public: 




string getName(); //{return 


"BC_NPI";} 


void init(); 




void initProbMax(); 




Void decay(EvtParticle *); 




private: 




EvtComplex Fpi(ql,q2); 




double pi3G(double Q2); 





Figure 3: The program structure. 

(modelName) is the name of the model that describes the specific decay, 
(pari) . . . (par m ) are the parameters of this model. The number of these parameters 
and their values depend on used model and specific decay. 

Modifying this .dec-file one can easily switch off the processes that are not currently needed 
and tune the processes that are used in the program. 

When Monte-Carlo generator is used, one needs to determine the maximal probability of 
the process. If this value is not set correctly, the efficiency of the program can be reduced, 
In EvtGen package maximum probability is determined by initProbMax method. If it is 
not set, the program would try to determine it automatically (simply generating 500 decays 
and taking the maximum value), and in the subsequent work any overflow incident would 
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be reported. 

The main function of any model is decay() method, where all physics of the process is 
coded. For a given decay kinematics this method calculates the amplitude of the process for 
all polarizations of initial and final particles. 

In order to include described in the previous section decays B c — > V(P) + nn the model 
BC_NPI was written. This model is realized in class EvtBcNPi and describes 5 c -meson 
decays both into vector and pseudoscalar mesons. The number of final 7r-mesons are limited 
to n < 3. The base class for both models is EvtDecayAmp (see figj3]), so polarizations 
of initial and final particles are taken into account accurately. The decayQ-methods are 
written according to relations (THj) . (TO)]) . (TTB]) from the previous section, and parameters of 
the models are read from determined by the user .dec-file. 

If we are describing B c — > V + mr decay, the corresponding record in the decay file should 
have the form 

(Br) (vec) (pi x ) . . . (pi n ) BC_NPI 
(maxProb) # maxProb 



where (Br) is the branching fraction of the decay, (vec) is the name of final vector particle, 
(pii) are the names of final 7r-mesons (pi+, pi- and piO for positively charged, negatively 
charge and neutral 7r-mesons respectively), and (Fi(0)), (c\) and (c|) are parameters of decay 
form factor Fi(q 2 ) (see eq.([5]) ). Note that lines in the decay file can be splitted freely and all 
symbols after # character are treated as comments. The order of the parameters, however, 
is important. In the case of the pseudoscalar meson in the final state the corresponding 
record in the decay file should have the form 

(Br) (vec) (ph) . . . (pi n ) BCLNPI 
(maxProb) # maxProb 



Note that in comparison with the previous example only one form-factor is left. 

To clarify this point let us consider the .dec file where only B c — > J/ipir + 7r° and B c — > 
B s n + are allowed. If SR form-factors set is used, this decay file should have the form 



(F^(0)> (cf) <cf > 
(Ff (0)> (c?+) <c^> 

(Mo)) {c\) {c\) 



# FAO 
#FA+ 
#FV 



(F+(0)> <cf+) (cD 



#F+ 
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Decay Bc+ 

# Be -> J/psi pi+ piO, SR form-factors set 

0.0017 J/psi pi+ piO BC.NPI 

330. # maxProp 

5.9 0.049 0.0015 # FAO 

-0.074 0.049 0.0015 # FAp 
0.11 0.049 0.0015; # FV 

# Be -> Bs pi+ , SR form-factors set 

0.18 B_sO pi+ BC.NPI 

250 #maxProb 
1.3 0.30 0.069; # Fp 

Enddecay 

IV. INSTALLATION INSTRUCTIONS AND COMPARISON WITH 

AVAILABLE DATA 

The EvtGen package requires CERNLIB, ROOT and CLHEP libraries to be properly 
installed. The location of the corresponding files can be altered by modifying the Makefile. 
Unfortunately the original package EvtGen is not supported now. Each experiment involved 
in the heavy quark research use its own branch of EvtGen for the decay generation. This 
why our module should be adapted to the concrete EvtGen branch. 

Here we present the variant of BC_NPI module which works with the last original version 



version of EvtGen 15]. This module version can be send on request by e-mail. 

The BC_NPI module is distributed as a source files EvtBcNpi.cc, EvtBcNpi.hh, sam- 
ple decay file be. dec, a simple test program test.ee, and a Makefile needed to include the 
considered model into EvtGen package. For installation of the model one should execute 
commands 

$ make lib 

$ sudo make install 

$ make clean 

Executing commands 
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$ make test 
$ ./test.exe 

one can run simple test program. This program generates 10 000 decays B c — > J/ip7i°7i + 
and saves a distribution over the invariant mass of the 7r°7r + -pair into file hist. root. The 
resulting historgam can be viewed with usual ROOT commands 

$ root hist. root 
hist->Draw() 

and should look like the distribution shown in figJU For comparison in this figure the results 
presented in [l3| are also shown. To remove the model from EvtGen package one can execute 
the command 

$ sudo make uninstall 

In addition, one should modify the decay-file and set the form-factor parameters. The 
format of this file was discussed in the previous section and parameters of B c — ^ J/ip and 
B c — > are presented in table fl] 

In Fig. H] we show the distribution over the invariant mass of light meson system in 
B c —7- J/ip + 2n decay is given in comparison with the prediction of work |l3j. One can 
easily see, that the results of Monte-Carlo are in good agreement with available data. For 
illustration we also show the same distribution with additional cut "p\ > 5 GeV 2 , where px is 
the transverse momentum of J/^-meson in the rest frame of initial B c . Such distributions, 
obviously, cannot be expressed in terms of spectral functions only. 



V. CONCLUSION 



_B c -mesons, that is particles build from b- and c-quarks are extremely interesting both from 
theoretical and experimental points of view. Since the masses of the constituent quarks are 
large in comparison with Aqcd, these particles give the opportunity to study QCD both in 
perturbative and nonperturbative regimes. Decays of I? c -mesons, on the other hand, can be 
caused only by weak interaction, so these particles give also information about electroweak 
sector of the Standard Model. Theoretically S c -mesons were studied rather thoroughly 
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Br 10" 




q 2 , GeV 2 



Figure 4: The distribution over the squared invariant mass q 2 



m* n in B c 



J/ip + 2ir decay 



generated within BC_NPI-model (histogram) in comparison with the predictions of 



131 ] (curve) . 



already, predictions for their masses, lifetimes and branching fractions of some, mainly two- 



body, decay modes are available. In our recent articles [13|, [14J we study also exclusive decays 
into a heavy quarkonia and a system of light mesons. In these articles the spectral function 
formalism is applied, when the phase space of this system is integrated out. It is clear, that 
main information about the dynamics of the processes is lost in this integration, so only 
branching fractions of the decays and distributions over the invariant mass of 7r-mesons can 
be obtained in the framework of this method. Additional information is required in the 
conditions of real experiments. 

Experimentally rather little is known currently about B c mesons. Only the mass of the 
ground state, its lifetime and branching fraction of two decay modes (B c — > J/ipiv and 
B c — > J/ipTr) were measured. According to theoretical estimates, a large yeld of B c - mesons 
is expected at LHC, so observation of other decay modes is possible. It is clear that, for the 



analysis of experimental data, it is necessary to have reliable theoretical predictions. 
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In this paper we consider the exclusive decays B c — > V(P)W* — > V(P) + mr, where 



1 < n < 3. Unlike our previous works 
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14|, we do not use the formalism of spectral 



functions, and describe in detail the process of transition of virtual IF-boson in a system of 
light mesons. This is realized with the help of the model for EvtGen generation package, that 
allows one to model these decays with different parametrizations of B c — > V(P) transition 
form-factors. This model is discussed in details in the present paper. 

This research is partially supported by Russian Foundation for Basic Research (grant 10- 
02-0006 la). The work of A. V. Luchinsky was also supported by non- commercial foundation 
"Dynasty" and the grant of the president of Russian Federation for young scientists with 
PhD degree (grant MK-406.2010.2). 
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